#include <iostream>
#include <fstream>
#include <cstdlib>
#include <cmath>
#include <ctime>

using namespace std;

class Two_Dim_Ising{

    public:

        Two_Dim_Ising(){
            srand(time(NULL));


            displayStr("Please Enter exchange coefficient::");
            cin >> J;

            displayStr("Please Enter the particle number in a chain:");
            cin >> chain_len;
            N = chain_len * chain_len;

            energy = -2 * N;
            energyAccumulation = 0;
            energySquareAccumulation = 0;

            magnetization = N;
            magnetizationAccumulation = 0;
            magnetizationSquareAccumulation = 0;

            displayStr("Please Enter the temperature of the system:");
            cin >> temperature;
            beta = 1.0 / temperature;

            displayStr("Please Enter the Running times:");
            cin >> times;

            lattice = new int*[chain_len];

            for (int i = 0; i < chain_len; i++){
                lattice[i] = new int[chain_len];
            }

            for (int x = 0; x < chain_len; x++){
                for (int y = 0; y < chain_len; y++){
                    lattice[x][y] = 1;
                }
            }

        }

        //Overload the constructor of  Two_Dim_Ising
        Two_Dim_Ising(double _B, double _J, int _chain_len, double t, long _times){
            srand(time(NULL));

			B = _B;
			J = _J;

            chain_len = _chain_len;
            N = chain_len * chain_len;

            //energy = -2 * N;
            energyAccumulation = 0;
            energySquareAccumulation = 0;

            //magnetization = N;
            magnetizationAccumulation = 0;
            magnetizationSquareAccumulation = 0;

            temperature = t;
            beta = 1.0 / temperature;

            times = _times;

            lattice = new int*[chain_len];

            for (int i = 0; i < chain_len; i++){
                lattice[i] = new int[chain_len];
            }

			acceptMove = 0;
			initializeToUp();

        }



        ~Two_Dim_Ising(){
            for (int i = 0; i < chain_len; i++){
                delete []lattice[i];
            }
            delete []lattice;
            lattice = NULL;
        }

        int run(){
            for (long i = 0; i < times; i++){
                step();
            }
        }

        double getSpin(int x, int y){
            if (x < 0 || x > chain_len || y < 0 || y > chain_len){
                cout << x << endl;
                cout << y << endl;
                displayStr("Out of Range!\n");
                exit(1);
            }
            return 1.0 * lattice[x][y];
        }

        double flipSpin(int x, int y){
            if (x < 0 || x > chain_len || y < 0 || y > chain_len){
                displayStr("Out of Range!\n");
                exit(1);
            }
            lattice[x][y] = -1 * lattice[x][y];
            return lattice[x][y];
        }

        double getEnergy(){
            return energyAccumulation / (times * N);
        }

        double getEnergySquare(){
            return energySquareAccumulation / (times * N);
        }

        double getHeatCapacity(){
            return (getEnergySquare() - getEnergy() * getEnergy()) / temperature;
        }

        double getMagnetization(){
            return magnetizationAccumulation / (times * N);
        }

        double getMagnetizationSquare(){
            return magnetizationSquareAccumulation / (times * N);
        }

        double getMagneticSusceptibility(){
            return (getMagnetizationSquare() - getMagnetization() * getMagnetization()) / temperature;
        }

        double getAccept(){
            return 1.0 * acceptMove / (times * N);
        }




    private:

		// exchange coefficient
		double J;

		// external magnetic field
		double B;

        // chain_len
        int chain_len;

        // number of particles
        int N;

        // system temperature
        double temperature;
        double beta;

        // system energy
        double energy;
        // system energy Accumulation
        double energyAccumulation;
        // system magnetization square Accumulation
        double energySquareAccumulation;


        //system magnetization
        double magnetization;
        // system magnetization Accumulation
        double magnetizationAccumulation;
        // system magnetization square Accumulation
        double magnetizationSquareAccumulation;




        //move accept
        long acceptMove;

        //lattice spin
        int **lattice;

        //run times
        long times;

        void displayStr(string str){
            cout << str;
        }

        // actually this function should not be defined here.I think that I need to write a class for generating random number
        double rand1(){
            double r = rand();
            if (r == RAND_MAX)
                return 0.0;
            else
                return 1.0 * r / RAND_MAX;
        }



		void initializeToUp(){
			magnetization = 0;
			for (int x = 0; x < chain_len; x++){
                for (int y = 0; y < chain_len; y++){
                    lattice[x][y] = 1;
					magnetization += lattice[x][y];
                }
            }
			getEnergyByLattice();
		}

		void getEnergyByLattice(){
			energy = 0;
			for (int x = 0; x < chain_len; x++){
				for (int y = 0; y < chain_len; y++){
					int down = y + 1;
					int right = x + 1;
					if (down  == chain_len) down = 0;
					if (right == chain_len) right = 0;
					energy += - B * lattice[x][y] + J *lattice[x][y] * (lattice[right][y] + lattice[x][down]);
				}
			}
			cout << energy << endl;
			cout << magnetization << endl;
		}

        void step(){
            for (int i = 0; i < N; i++){

                int x = (int)(1.0 * rand1() * chain_len);
                int y = (int)(1.0 * rand1() * chain_len);

                int xNext = x + 1;
                int xBefore = x - 1;

                int yNext = y + 1;
                int yBefore = y - 1;

                if (xNext == chain_len) xNext = 0;
                if (yNext == chain_len) yNext = 0;

                if (xBefore == -1) xBefore = chain_len - 1;
                if (yBefore == -1) yBefore = chain_len - 1;

                double energyBefore = J  * getSpin(x, y) * (getSpin(xNext, y) + getSpin(xBefore, y) + getSpin(x, yBefore) + getSpin(x, yNext)) - B * lattice[x][y];
                double energyAfter = -1.0 * J * getSpin(x, y) * (getSpin(xNext, y) + getSpin(xBefore, y) + getSpin(x, yBefore) + getSpin(x, yNext)) + B * lattice[x][y];
                double de = energyAfter - energyBefore;

                if (de <= 0 || exp(-1 * beta * de) > rand1()){
                    double dm = flipSpin(x, y);
                    acceptMove++;
                    energy += de;
                    magnetization += 2 * dm;
                }

                energyAccumulation += energy;
                energySquareAccumulation += energy * energy;
                magnetizationAccumulation += magnetization;
                magnetizationSquareAccumulation += magnetization * magnetization;


            }
        }

};


int main(){


    ofstream ostream;
    ostream.open("two_dim_ising_model.txt");


    for (double t = 0.1; t < 10.0; t += 0.1){
        Two_Dim_Ising *x = new Two_Dim_Ising(0, -1, 16, t, 10000);
        x->run();
        ostream << t << "\t" << x->getEnergy() << "\t" << x->getHeatCapacity() << "\t" << x->getMagnetization() << "\t" << x->getMagneticSusceptibility() << "\t" << x->getAccept() << endl;
        delete x;

    }

    ostream.close();
    return 0;
}
